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Abstract. We review recent progress in modeling the probability distribution of wave 
heights in the deep ocean as a function of a small number of parameters describing 
the local sea state. Both linear and nonlinear mechanisms of rogue wave formation are 
considered. First, we show that when the average wave steepness is small and nonlinear 
wave effects are subleading, the wave height distribution is well explained by a single 
"freak index" parameter, which describes the strength of (linear) wave scattering by 
random currents relative to the angular spread of the incoming random sea. When the 
average steepness is large, the wave height distribution takes a very similar functional 
form, but the key variables determining the probability distribution are the steepness, 
and the angular and frequency spread of the incoming waves. Finally, even greater 
probability of extreme wave formation is predicted when linear and nonlinear effects 
are acting together. 



Linear and Nonlinear Rogue Wave Statistics in the Presence of Random Currents 2 
1. Introduction 

Tales of freak waves by lucky survivors used to be taken with a large grain of salt. Were 
the sailors making excuses for bad seamanship? The first such wave to be measured 
directly was the famous New Year's wave in 1995 [1]. With modern cameras and video, 
not to mention satellites [21 13] , it is no longer controversial that freak or rogue extreme 
waves exist on the world's great oceans [U |5l |6]. 

Any realistic seaway (an irregular, moderate to rough sea) is comprised of a 
superposition of waves differing in wavelength and direction, with random relative 
phases. Supposing that the dispersion in wavelengths is not too large, and assuming 
uniform sampling, Longuet-Higgins ^ exploited the central limit theorem to derive a 
large number of statistical properties of such wave superpositions, including of course 
wave height distributions. From this viewpoint, extreme waves are the result of unlucky 
coherent addition of plane waves corresponding to the tail of the Gaussian distribution 
(see Eq. (|9]) below). As explained below, wave heights greater than about 4 a in the tail 
of the Gaussian are classified as extreme. The problem has become how to explain why 
the observed number of rogue wave events is greater than the number A a out in the 
Longuet-Higgins theory. 

For the following discussion, it is important to understand why a 20 meter wave in a 
sea where the significant wave height (SWH, defined as the average height of the highest 
one third of the waves) is 18 meters is far less onerous than a 20 meter wave where the 
SWH is 8 meters. An established seaway of uniform energy density (uniform if averaged 
over an area large compared to the typical wavelength) is "accommodated" over time 
and distance, through nonlinear energy transfer mechanisms. Seaways of higher energy 
density develop correspondingly longer wave periods and wavelengths, even with no 
further wind forcing, keeping wave steepness under control as a result. 

This accommodation process is one of the ways nonlinear processes are implicitly 
lurking behind "linear" theories, in that the input into the linear theories, i.e., the 
SWH, the period, dispersion in direction, and dispersion in wavelength are all the result 
of prior nonlinear processes. A 20 meter wave in a sea of SWH 8 meters is necessarily 
very steep, possibly breaking, with a deep narrow trough preceding it. The tendency 
for steep waves to break is an often devastating blow just as the ship is sailing over an 
unusually deep trough before meeting the crest. 

Observational evidence has shown that the linear Longuet-Higgins theory is 
too simplistic Recent advances in technology have allowed multiple wave tank 
experiments and field observations to be conducted, confirming the need for a more 
realistic theory to explain the results [HI |9]. An obvious correction is to incorporate 
nonlinear wave evolution at every stage, rather than split the process into an implicit 
nonlinear preparation of the seaway followed by linear propagation. Clearly the exact 
evolution is always nonlinear to some extent, but the key is to introduce nonlinearities 
at the right moment and in an insightful and computable way. Realistic fully nonlinear 
computations wave by wave over large areas are very challenging, but initial attempts 
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have been made to simulate the ocean surface using the full Euler equation both on 
large scales |T0] and over smaller areas |Tll [12] . 

Surprisingly, investigation of nonlinear effects is actually not the next logical step 
needed to improve upon the Longuet-Higgins model. Indeed the full linear statistical 
theory had not been given, for the reason that uniform sampling assumed by Longuet- 
Higgins is not justified. A nonuniform sampling theory, which does not assume that the 
energy density is uniformly distributed over all space, is possible and is still "linear." 
Moreover, the parameters governing a nonuniform sampling are knowable. Inspired by 
the work of White and Fornberg [13] , the present authors showed that current eddies 
commonly present in the oceans are sufficient to cause the time-averaged wave intensity 
to be spatially non-uniform, and to exhibit "hot spots" and "cold spots" some tens 
to hundreds of kilometers down flow from the eddies. We emphasize that in terms 
of wave evolution, the refraction leading to the patchy energy density is purely linear 
evolution. The key ideas are (1) that waves suddenly entering a high energy patch 
are not accommodated to it and grow steep and dangerous, and (2) the process is still 
probabilistic and the central limit theorem still applies, with the appropriate sampling 
over a nonuniform distribution. The high-energy patches will skew the tails of the 
wave height distribution, perhaps by orders of magnitude. This was the main point in 
reference |14j . 

There is no denying the importance of nonlinear effects in wave evolution, and a full 
theory should certainly include them. On the other hand a nonlinear theory that fails to 
account for patchy energy density is missing an important, even crucial effect. The linear 
theory needs to be supplemented by nonlinear effects however, since the accommodation 
of the waves to the presence of patchy energy density needs to be considered. It is our 
goal here to review progress along these lines and point the way to a more complete 
theory. We first review the nonuniform sampling linear theory and then discuss newer 
simulations using the nonlinear Schrodinger equation (NLSE). Finally, we show that 
even larger rogue wave formation probabilities are predicted when linear and nonlinear 
formation mechanisms are acting in concert. 

Rogue wave modeling would benefit greatly from a comprehensive, accurate, and 
unbiased global record of extreme wave events, supplemented by data on local ocean 
conditions, including current strength, SWH, steepness, and the angular and spectral 
spread of the sea state. Such a record, not available at present, would allow for direct 
statistical tests of linear and nonlinear theories of rogue wave formation. Anecdotal 
evidence does suggest that rogue waves may be especially prevalent in regions of strong 
current, including the Gulf Stream, the Kuroshio Current, and especially the Agulhas 
Current off the coast of South Africa. Consequently, the Agulhas Current in particular 
has attracted much attention in rogue wave research [5l[T5]. However, anecdotal evidence 
from ships and even oil platform measurements cannot provide a systematic, unbiased, 
and statistically valid record that would support a correlation between possibly relevant 
variables and rogue wave formation probability. Instead, satellite-based synthetic 
aperture radar (SAR) [21 [3] is currently the only method that shows potential for 



Linear and Nonlinear Rogue Wave Statistics in the Presence of Random Currents 4 



monitoring the ocean globally with single-wave resolution, but validating the surface 
elevations obtained by SAR is a challenge. The SAR imaging mechanism is nonlinear, 
and may yield a distorted image of the ocean wave field; the nonlinearity is of course 
of particular concern for extreme events [16] . Recently, an empirical approach has been 
proposed that may accurately obtain parameters such as the SWH from SAR data [17] . 

2. Linear Wave Model in Presence of Currents 

2.1. Ray Density Statistics 

To understand the physics of linear rogue wave formation in the presence of currents, it 
is very helpful to begin with a ray, or eikonal, approximation for wave evolution in the 
ocean [131 [H], 

dk du{f,k) df duj{f,k) . 

J[= ar ' ~dt^ dk ' 

where r is the ray position, k is the wave vector, and u is the frequency. For surface 
gravity waves in deep water, the dispersion relation is 

uj{r,k) = ^J^\ + k-U{f), (2) 

— * 

where U{r) is the current velocity, assumed for simplicity to be time-independent, and 
g = 9.81 m/s is the acceleration due to gravity. The validity of the ray approximation 
depends firstly on the condition \k\^ ^ 1, where ^ is the length scale on which the 
current field U{r) is varying, physically corresponding to the typical eddy size. This 
condition is well satisfied in nature, since wave numbers of interest in the deep ocean 
are normally of order k ~ 27r/(100m), while the typical eddy size may be ^ ~ 5 km or 
larger. Secondly, scattering of the waves by currents is assumed to be weak, i.e., the 
second term in equation ([2]) should be small compared to the free term. This again is 
well justified since eddy current speeds \U\ are normally less than 0.5 m/s, whereas the 
wave speed v = du/dk ^J gjAk is greater than 5 m/s. In section [23] below, we will 
explicitly compare the ray predictions with results obtained by exact integration of the 
corresponding wave equation. 

In the numerical simulations shown in figure [H we follow White and Fornberg [13] 
in considering a random incompressible current field in two dimensions, with zero mean 
current velocity, generated as 

Vx{?) = —dip{f)/dy; Uy{f) = dip{r)/dx (3) 

from the scalar stream function ilj{r). The stream function itself is Gaussian distributed 
with Gaussian decay of spatial correlations: 

W) = ; W)W) ~ e-(^-"")'/2«' , (4) 

and the overall current strength is described by u"^^^ = |f/(r)p. The specific choice 
of a Gaussian distribution for the stream function is made for convenience only. The 
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Figure 1. A ray density map I{x,y) is calculated for rays moving through a 640 km by 
640 km random eddy field, with rms eddy current Wrms — 0-5 m/s and eddy correlation 
length ^ = 20 km. Here bright regions represent high density. The rays are initially 
distributed uniformly along the left edge of each panel, with angular spread A9 around 
the +x (rightward) direction, and with frequency uj = 27r/(10sec), corresponding to 
velocity v — 7.81 m/s in the absence of currents. The left and right panels illustrate 
A9 = 10° and AO = 20°, respectively. 

detailed structure of the individual eddies on the scale ^ has no effect on the final rogue 
wave statistics as long as the current is weak ("Urms ^ v), since each ray must travel a 
distance ^ C, before being appreciably scattered. Each panel in figure [H represents a 
640 km by 640 km random eddy field, with rms eddy current Wrms = 0.5 m/s and eddy 
correlation length ^ = 20 km. 

The initial swell, entering from the left in each panel, is characterized by a single 
frequency u = 27r/(10sec) (and thus a single wave number k = u^/g = 0.04 m~^ and 
a single wave speed v = \J g/Ak = 7.81 m/s). As discussed in Ref. [11], within the 
context of a linear model, a nonzero frequency spread affects rogue wave formation only 
at second order in the spread Au, and may be neglected for all practical purposes. In 
contrast, the angular spread of the incoming sea is very important in determining rogue 
wave statistics. In this figure, we consider an initially Gaussian angular distribution 
p{9) ~ g-e^/aCAe)^^ where 9 is the wave vector direction relative to the mean direction 
of wave propagation. Here all rays begin at the left edge of each panel, uniformly 
distributed in the y direction, and the mean direction of wave propagation is rightward. 
The left and right panels illustrate the behavior for two different values of the initial 
angular spread A9. 

In both panels we observe bright streaks or branches, corresponding to regions of 
larger than average ray density I{x, y), and thus larger than average wave intensity. The 
branches may be understood by considering briefly the limiting (unphysical) case of a 
unidirectional initial sea state {A9 = 0), corresponding to a single incoming plane wave. 
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In the ray picture, and in the coordinates of figure [H the initial conditions are in this 
hmit characterized by a one-dimensional phase space manifold (x, y, k^, ky) = (0, y, k, 0), 
where k is the fixed wave number, and y varies over all space. As this incoming plane 
wave travels through the random current field, it undergoes small-angle scattering, with 
scattering angle ~ Urms/v after traveling one correlation length ^ in the forward direction. 
Eventually, singularities appear that are characterized in the surface of section map 
[l/(0), fcj^(O)] — 7- [y{x),ky{x)] by 5y{x) / 5y{Q) = 0, i.e., by local focusing of the manifold 
of initial conditions at a point {x,y). 

The currents leading to such a focusing singularity may be thought of as forming 
a 'bad lens.' Whereas a lens without aberration focuses all parallel incoming rays to 
one point, a bad lens only focuses at each point an infinitesimal neighborhood of nearby 
rays, so that different neighborhoods get focused at different places as the phase-space 
manifold evolves forward in x, resulting in lines, or branches, of singularities. The typical 
pattern is an isolated cusp singularity, S'^x{y)/Sx{Oy = 0, followed by two branches of 
fold singularities, as shown in figure El 




Figure 2. A cusp singularity, followed by two branches of fold singularities, is formed 
as initially parallel rays pass through a focusing region. The two branches appear 
because the focal distance varies with the distance of approach from the center, as in 
a 'bad' lens with strong spherical aberration. After averaging over incident directions, 
the singularities will be softened but not washed away completely [2] 

A simple scaling argument [131 [13 HH] shows that the first singularities occur after 
a median distance ?/ = L ~ ^{uj:ms/v)~'^^^ ^ ^ along the direction of travel, when the 
typical ray excursion in the transverse x direction becomes of order ^. Thus, each ray 
passes through many uncorrelated eddies before a singularity occurs, and a statistical 
description is well justified. For realistic parameters, L ~ 100 km or more is typical. 
The cusp singularities formed in this way are separated by a typical distance ^ in the 
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transverse direction, and thus the rms deflection angle by the time these singularities 
appear scales as 

50~e/L~K^sAO'^'- (5) 

We note that the typical deflection angle 69 does not depend on the eddy size but 
only on the velocity ratio u^^s/v: faster currents cause larger deflection. For the 
input parameters used in flgure [1], the median distance to the flrst singularity is 
L = 7.5C, = 150 km, and the rms deflection at the point of singularity is 66 = 18°. 




Figure 3. Predicted tsunami wave heights from the Tohoku earthquake, a 9.0 
magnitude undersea earthquake that occurred on March 11, 2011, off the coast of 
Japan. A branching structure is clearly visible as the waves move outward from the 
epicenter. (Source: NOAA Center for Tsunami Research.) 



Similar phenomenology can give rise to wave focusing and rogue wave formation in 
shallow water, where the dispersion relation of equation ([2]) is replaced with u{r, k) = 
gktanh{kh{r)), and varying depth h{f) takes the place of the varying current U{r) as 
the origin of scattering [20j. The same mechanism can lead to ampliflcation of tsunami 
waves [m [22] where because of the long wavelength, shallow water equations apply. 
Fig. in] shows a striking recent example of a predicted tsunami wave height map, in 
which the branched flow structure is unmistakably present. More generally, singularities 
and branched flow due to focusing in random media have been investigated in contexts 
as diverse as electron flow in a two-dimensional electron gas [23], ocean acoustics [21], 
twinkling of starlight [25], and rain shower activation in turbulent clouds [26]. Recently, 
universal expressions have been obtained describing the branching statistics for a large 
class of such systems, and valid at all distances from a source |l9j. 

For finite initial angular spread A6, the singularities are softened, and the finite 
contrast between the peak ray density in the branches and the background intensity is 
governed for ^ 1 and 56* ^ 1 by the ratio 

2/3 
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Figure 4. The ray density distribution, for an initial sea state of uniform density 
scattered by a random eddy current field, is shown for several values of the freak index 
7. The input parameters are chosen as in figure[l] with initial angular spread A9 — 25°, 
15°, and 5° corresponding to freak index 7 = 0.72, 1.20, and 3.60, respectively. The 
mean intensity is normalized to unity in each case. The dashed curves are fits to the 
distribution of Eq. ([7]). 



which we refer to as the freak index [m]. Of particular interest is the regime of small 7, 
where the scattering characterized by 66 is weak compared to the initial angular spread 
of the incoming sea. In this limit, the scattering produces only small perturbations 
of order 7"^ in the ray density I{x,y), in units where the initial (uniform) density is 
-^0 = 1 [E] ■ Specifically, as seen in figure HI the distribution of ray intensities in this 
regime may be well described by a distribution 



where the number of degrees of freedom scales with the freak index as 7"^. The 
proportionality constant may be obtained numerically by a fit to the data: 
« 45 

N=- = -. 8 

7^ 7^^ 

In the limit 7 — associated with zero current, we have — t- 00, and we recover as 
expected the uniform ray density distribution g{I) = S{I — 1). 

2.2. Implications for Wave Statistics 

In the Longuet-Higgins random seas model the sea surface elevation above the 
average elevation is given by Iie({x, y,t), where C is a random superposition of many 
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plane waves with differing directions and frequencies. By the central limit theorem, 
( is distributed as a complex Gaussian random variable with standard deviation a. 
Furthermore, for a narrow-banded spectrum {6co <C uj) the wave crest height H is equal 
to the wave function amplitude \(\, and the probability of encountering a wave crest of 
height H or larger is 

PRayleigh(i^) = . (9) 

Due to an exact symmetry between crests and troughs in a linear wave model, a crest 
height of H corresponds to a wave height (crest to trough) of 2H. Conventionally, a 
rogue wave is defined as 2H > 2.2 SWH, where the significant wave height SWH is 
the average of the largest one third of wave heights in a time series, or approximately 
SWH ^ 4.0(7. Thus the condition for a rogue wave is if > 4Aa, and the random 
seas model predicts such waves to occur with probability -PRayieigh(4.4cT) = 6.3 ■ 10~^. 
Similarly, extreme rogue waves may be defined by the condition 2H > 3 SWH or 
H > 6. Oct, and these are predicted to occur with probability -PRayieigh(6.0(T) = 1.5 • 10~^ 
within the random seas model. As discussed in section [H the random seas model 
greatly underestimates the occurrence probability of extreme waves, when compared 
with observational data [2]. 

What are the implications of scattering by currents, as discussed in section 12. 
on the wave height statistics? Within the regime of validity of the ray approximation, 
we have at any spatial point (x, y) correspondence between the ray density I{x, y) and 
the wave intensity = \({x,y,t)\'^, averaged over time. Thus, in contrast with the 
original Longuet-Higgins model, the time-averaged wave intensity is not uniform over 
all space but instead exhibits "hot spots" and "cold spots" associated with focusing 
and defocusing in the corresponding ray equations. At each point in space (assuming of 
course that the currents are stationary), the central limit theorem and thus the Rayleigh 
distribution still apply, and we have 

P(.,,)(i7) = e-^'/2'^'^(^'^), (10) 

where I{x, y) is the local ray density, normalized so that the spatial average is unity, and 
cr^ is the variance of the surface elevation in the incoming sea state, before scattering 
by currents. This is the situation a ship experiences at a given position. 

Now averaging over space, or over an ensemble of random eddy fields with a given 
rms current speed, we obtain a total cumulative height distribution 

PtotdH)= / rf/^?(/)e-^'/2'^'^ (11) 

In equation (ITTl) . the full cumulative distribution of wave heights for a given sea state has 
been expressed as a convolution of two factors: (i) the local density distribution g{I), 
which can be extracted from the ray dynamics, and (ii) the universal Longuet-Higgins 
distribution of wave heights for a given local density. Similar decompositions of chaotic 
wave function statistics into non-universal and universal components have found broad 
applicability in quantum chaos, including for example in the theory of scars [281 EH]. 
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the context of rogue waves, a similar approach was adopted by Regev et al. to study wave 
statistics in a one-dimensional inhomogeneous sea, where the inhomogeneity arises from 
the interaction of an initially homogeneous sea with a (deterministic) long swell [30] . 

Using the previously obtained ray density distribution in the presence of currents, 
equation ([7]), we obtain the K-distribution [31] 



JV 
2 



Ptotai(^) = 2-^;7i^777^v^i^7v/2 (VNHo) , (12) 



JV 
2 



Ptotai(x) = 2^,^^^j^Kn/2 { 'iVNx ) , (13) 



T{N/2) ""'n J 

where Kn{y) is a modified Bessel function. 

Defining the dimensionless variable x = 2if/SWH ^ 2H/ (4cr), so that a rogue wave 
is given by x = 2.2 and an extreme rogue wave by x = 3.0, we find the probability of a 
wave height exceeding x significant wave heights: 

(^y/Nxy 

T{N/2) 

to be compared with the random seas prediction 

i'Rayleigh(a;) = C'^^' (14) 

in the same dimensionless units. We recall that N in equation (1121) or (fT3|) is a function 
of the freak index 7, as given by equation (IHl). 

To examine the predicted enhancement in the probability of rogue wave formation, 
as compared with random seas model ([9]), we may consider two limiting cases. Keeping 
the wave height of interest fixed, and taking the limit 7 — )■ 0, i.e. N — )• 00, we obtain 
the perturbative result 

4 



-Pperturb \X) 



l + -{x'-x') 

l + ^{x'-x') 
b 



-PRayleigh(a;) (15) 
-PRayleigh(a;) , (16) 



valid for x^ <^ N, or equivalently x^7 ^ 1. Thus, in the limit of small freak index, the 
distribution reduces, as expected, to the prediction of the random seas model. Analogous 
perturbative corrections appear for quantum wave function intensity distributions in the 
presence of weak disorder or weak scarring by periodic orbits [321 [33] . 

Much more dramatic enhancement is observed if we consider the tail of the intensity 
distribution (x — >■ oo) for a given set of sea conditions (fixed 7 or A^). Then for x ^ N^^'^, 
or equivalently x^^ ^ 1, we obtain the asymptotic form 

V-l 
2 



Nx 

-^asymptotic f-^) ^TTTTTTTT 6 



r(A^/2) 



Nx 



JV-1 
2 
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i.e., the probability enhancement over the random seas model is manifestly super- 
exponential in the wave height x. 

Predicted enhancements in the probability of rogue wave and extreme rogue wave 
formation, based on equations (fTSll and ([HD, are shown in table [H We notice in 
particular that enhancements of an order of magnitude or more are predicted in the 
extreme tail, even for moderate values of the input parameters, corresponding to 
7 ~ 0.72 - 1.2 or ~ 30 - 85. 





7 
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22 
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2.0 


9.8 
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0.60 


125 


1.7 


5.7 



Table 1. The N parameter of equation ([T3|) . and the associated enhancement in 
the probability of rogue wave formation (wave height 2H = 2.2 SWH) as well as 
the enhancement of the probability of extreme rogue wave formation (wave height 
2H — 3.0 SWH) are calculated for several values of the incoming angular spread A9 
using equations ([6|), (|8]), and ([T3l) . Here E{x) — /totai(a;)/^R,ayicigh(a;). In all cases 
we fix the rms current speed Urms — 0.5 m/s and mean wave speed v ~ 7.8 m/s, so 
^6* = 18°. 



2.3. Numerical Results for Linear Wave Equation 



The theoretical predictions of equation (]T3l) are based on several approximations, 
including the assumption of local Rayleigh statistics. To see whether the approximations 
we have made are valid, we compare the theoretical predictions with direct numerical 
integration of the current-modified linear Schrodinger equation, which is obtained from 
the third-order current-modified nonlinear Schrodinger (CNLS) equation [M] by setting 
the nonlinear term to zero. CNLS governs the modulations of weakly nonlinear water 
waves around a mean frequency and mean wave vector, incorporating the effect of 
currents, and is presented in full in Sec. 13. II below. In dimensionless variables, the linear 
equation for the wave envelope describing the wave modulations is expressed as [31] 



iAr - Iaxx + -Ayy - hU^A = . (18) 

o 4 

Here A{X, Y, T) is the wave envelope, defined by separating out the carrier wave 
propagating with mean wave vector k = kox, 

C{X, Y, T) = koA{X, Y, T)e''"'^-'^' 

= koA{X,Y,T)e'''-'^/\ (19) 
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and 

(X, Y, T) = (kox - ^ V^t, koy, x/^ot) (20) 

are dimensionless space and time coordinates. We also note that Eq. fllSp may be 
obtained directly from the dispersion relation ([2]), by expanding u and k around cuq and 
kox, respectively. 
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Figure 5. The probability of exceeding wave height 2H, in units of the significant 
wave height SWH, is shown for an incoming wave speed v — 7.8 m/s, incoming angular 
spread = 5.7°, and rms current speed Umis — 0.5 m/s. The solid curve shows the 
results of a numerical simulation performed on a 20 km by 40 km field, with typical 
eddy size f — 800 m, while the dashed curve represents equation ([T3| with N ~ 6.8. 
The Rayleigh (random seas) prediction of equation (|14p is shown for comparison. 



The calculations are performed on a rectangular field measuring 40 km along the 
mean direction of propagation and 20 km in the transverse direction, with typical eddy 
size ^ = 800 m. (We note that a very small value for the eddy size is chosen to 
maximize the statistics collected; this is also a "worst case" scenario for the theory, as 
the ray approximation is expected to work ever better as the ratio of eddy size to the 
wavelength increases.) Equation ( fT8l) is integrated numerically using a split-operator 
Fourier transform method [35], on a 1024 by 512 grid. The incoming wave is a random 
superposition of a large number of monochromatic waves with directions uniformly 
distributed around the mean direction ^ = with standard deviation A6. Without loss 
of generality, the incoming wave number is fixed at ko = 27r/(156m), corresponding to 
a frequency ojq = y/gko = 27r/(10sec) and a group velocity v = 7.81 m/s. Each run 
simulates wave evolution for 5 • 10^ sec or 5 • lO'' wave periods, sufficient for the wave 
height statistics to converge. 

The results for A6 = 5.7°, corresponding to a very large freak index 7 = 3.15, 
are shown in figure [5l The results are compared both with the theoretical prediction 
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of equation (1131) (here N = 6.8) and with the basehne Rayleigh distribution of 
equation f|T^ . This is an extreme scenario, in which the occurrence probabihty of 
extreme rogue waves (3 times the significant wave height) is enhanced by more than three 
orders of magnitude. Even better agreement with the theoretical model of equation (fT3|) 
obtains for more moderate values of 7, corresponding to larger A^. 
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Figure 6. The wave height distribution for a random incoming sea scattered by 
random currents is obtained numerically for four values of the rms current speed 
Mims and four values of the incoming angular spread A9. In each case, a fit to the 
K-distribution (equation (1121) ) yields the number-of-degrees-of- freedom parameter TV, 
(describing deviations from Rayleigh statistics), which is plotted as a function of the 
freak index 7 (defined in equation As in the previous figures, the wave speed is 
fixed a,t V = 7.81 m/s. The solid line is the theoretical prediction of equation ([8)). 



In figure [6] we repeat the numerical simulation for four different values of the 
incoming angular spread and four different values of the rms current speed Urms- 
In each case, the numerically obtained wave height distribution is fit to a K-distribution 
(equation ( |T2l) ). and the resulting value of (which fully describes the strength of 
deviations from Longuet-Higgins statistics) is plotted as a function of the freak index 
7. Excellent agreement is observed with the power-law prediction of equation ([8]) all 
the way up to 7 ~ 2 (corresponding to ~ 10), even though the analytic prediction 
was obtained in a small-7 approximation. The regime in which the analytic formula 
(IH]) works well includes most conditions likely to be found in nature (e.g., all but the 
first row of tabled]). Referring again to tabled! we observe that the theory accurately 
describes enhancements of up to three orders of magnitude in the formation probability 
of extreme rogue waves. Modest deviations from the analytic formula are observed 
numerically at very large values of 7 (corresponding to even larger enhancements). 
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Figure 7. The left panel shows the distribution of the experimentally measured time- 
averaged microwave intensity s, found for 780 positions in the random wave field, is 
compared with the distribution of equation ([7]), with N = 32. The inset shows the 
same data on a logarithmic scale. The right panel shows a time series of the wave 
intensity at a single point shows a "rogue wave" event. The inset shows a snapshot of 
the wave intensity field near this point at the moment of the extreme event. 

2.4- Experimental Demonstration of Linear Rogue Wave Formation 

Direct experimental verification in the ocean of tlie statistical predictions made 
analytically in section 12.21 and confirmed numerically in section 12.31 is obviously highly 
desirable. Unfortunately no observational data set exists at this time that would allow 
the tail of the wave height distribution to be studied as a function of the freak index 
7, i.e., as a function of the rms current speed Wrms and of the angular spread A^. 
Recently, however, experiments in open quasi-two-dimensional microwave cavities with 
disorder [2^ have found a strong enhancement in the occurrence probability of high- 
amplitude waves, which may be interpreted as "rogue waves" in this analog system. 
In the microwave system, randomly placed brass cones play the role of random ocean 
currents, and a movable source antenna enables incoming waves to arrive from different 
directions. A movable drain antenna acts a weak probe, and allows for a spatial 
mapping of the wave fields within the scattering arrangement. A great advantage of 
the microwave system is that the electromagnetic wave equation is linear, so that the 
observed enhancement in the tail of the wave height distribution may serve in principle 
as a direct experimental test of the theory developed in the previous sections. 

In the left panel of figure [TJ the time-averaged wave intensity s is found for different 
positions of the probe, and the probabihty distribution g{s) is shown. Most of the 
distribution is well described by the distribution of equation (I7l). We note that in the 
absence of disorder, the time-averaged intensity would be position-independent and g{s) 
would reduce to 6{s — 1) (A^ — )■ oo in equation ([7])). The inset in the left panel shows 
additional rare events in the far tail, that are not described by the distribution |27j . 
The right panel in figure [7] shows time series data of the wave intensity at a single 
point, including an extreme event observed in the experiment, and the inset shows a 
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snapshot of the wave intensity in the region at the moment corresponding to this extreme 
event. The event presented here has wave height 2H = 5.3 SWH, and events of this 
magnitude of greater are observed with probabihty 1.3 x 10"^ in the experiment, which 
is an enhancement of 15 orders of magnitude compared to the Rayleigh distribution. 

These results confirm that hnear scattering is a sufficient mechanism for a large 
enhancement in the tail of the wave height distribution, even when nonlinearity is 
entirely absent from the physical system being studied. 

3. Nonlinear Wave Model 

We have already seen (e.g., in table [I]) that under physically realistic sea conditions, 
linear wave dynamics, with nonlinearity only in the corresponding ray equations, 
are sufficient to enhance the incidence of extreme rogue waves by several orders of 
magnitude. At the same time, the true equations for ocean wave evolution are certainly 
nonlinear, and furthermore the nonlinear terms, which scale as powers of the wave height, 
manifestly become ever more important in the tail of the wave height distribution. Thus, 
a fully quantitative theory of rogue wave statistics must necessarily include nonlinear 
effects, which we address in the following. 

3.1. Nonlinear Schrddinger Equation 

The original Nonlinear Schrodinger Equation (NLSE) for surface gravity waves in deep 
water was derived by Zakharov using a spectral method [36], and is valid to third 
order in the steepness e = k^H, where H is the mean wave height. Subsequently, 
the NLSE was extended to fourth order in e by Dysthe [32| and then to higher order 
in the bandwidth Au/u by Trulsen and Dysthe [38] • The Trulsen-Dysthe equations 
include frequency downshifting [39], the experimentally observed reduction in average 
frequency over time [ID]; however the physics of frequency downshifting may not yet be 
fully understood [H]. 

In our simulations we implement the current-modified 0{e'^) NLSE, as derived by 
Stocker and Peregrine in dimensionless form |34j : 



where the the linear and third-order terms are collected on the left hand side of 
equation fl2T]) . Here $, $c; and B represent the mean flow, surface current, and 
oscillatory parts, respectively of the velocity potential 0: 



where the second-harmonic term B2 is function of B and its derivatives, {X, Y, T) 
are dimensionless space and time coordinates defined previously in equation (!20|) . and 



ij^^ — (^Bxx — '^Byy) — 7:B\B\^ — B^cX = ypiBxxx — QByyx) 

o Z Lb 

+ ^xB + '-B{BB*x - QB*Bx) + ti^^cXT - ^cz)B - iWn^, ■ WhB (21) 




(22) 
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6 = k^x — y/gkot = X — T/2 is the phase. The surface elevation, which is the quantity 
of interest for our purposes, is similarly expanded as 



C = h-' 



C + Cc + l {A'' + A2e^'' + A^e'"' + c.c.) 



(23) 



where the expansion coefficients may be obtained from the velocity potential as 

A =iB + -^B^ + 4o(5xx - 25„^) + 

2kQ %kl ^ '''' ra; ^ g I I 

A2= - -B^ + ^BB, (24) 
2 ko 



3i 



?3 



= B' 

^ 8 

***Here both B and A are of order e, and proportion to e. By changing the magnitude 
of B or A in the incoming wave, we can set steepness to different value. 

In the simulatin, for the simplicity, we works in the frame of reference moving 
with the velocity vq = (co + Uo,Vo), so and $c in equation fl22|) is zero. The 
incoming wave is a random superposition of a large number of monochromatic waves 
with different frequencies and propagating directions. Thus the initial wave could be 
prepared analytically, as a linear summation of a large number of plain wave, 

N N 

^(f,t) = ^0, = ^A,e*-^^ (25) 
1 1 

where ki is the random wave vector for each monochromatic wave. For our setup, the 
wave vector can be expressed as 

k = {ko + k') ■ {cose' X + sine' y) (26) 

where k' is a random variation in wave number follows a Gaussian distribution whose 
half height width is Afc, and e' is the angular spread which is a normal distribution with 
stand deviation A6'. 

In the following examples, equation ( 12T1) is integrated numerically with the current 
set to zero, in order to investigate systematically and quantitatively the effect of 
nonlinear focusing. In nature, the interplay between linear and nonlinear mechanisms is 
also of great interest, and may give rise to even stronger enhancement in the probability 
of rogue wave occurrence than either effect individually, as demonstrated below in 
section H (see also [i2lli8]). 



3.2. Height Distribution 

As in the linear case, the split-operator Fourier transform method is used to integrate 
equation fl2T]) numerically. The rectangular field measuring 20 km along the mean 
direction of propagation and 10 km in the transverse direction is discretized using a 1024 
by 512 grid. The incoming state is a random superposition of plane waves with wave 
numbers normally distributed around k^ with standard deviation A/c, and directions 
uniformly distributed around the mean direction ^ = with standard deviation A^. 
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Without loss of generality we fix the mean incoming wave number at kg = 27r/(156m), 
as in section [2l The steepness koH is adjusted by varying the mean height H of the 
incoming sea. Each run simulates wave evolution for 4 ■ 10® sec or 4 ■ 10^ wave periods. 

Typical results are represented by solid curves in figure [H where we fix Ak/ko = 0.1 
and A6 = 2.6° (as we will see below, the values of A6 required to see very strong 
effects from nonlinear focusing are typically smaller than those needed to observe 
significant deviations from Rayleigh by linear scattering). The cumulative probability 
distribution of the wave height 2H, in units of the significant wave height SWH, is 
shown for four nonzero values of the wave steepness e. As expected, the Rayleigh 
probability distribution of equation fll4p is recovered in the limit e — )■ 0, and ever stronger 
enhancement in the tail is observed as the steepness of the incoming sea increases. The 
occurrence probability of extreme rogue waves, 2if/SWH = 3.0, is enhanced by one to 
three orders of magnitude for the parameters shown. 

To understand the functional form of the distributions in figure [HI we again make 
use of the local Rayleigh approximation discussed above in section 12.21 Here the wave 
height distribution is given locally in space and time by a Rayleigh distribution around 
the local mean height (corresponding to a locally random superposition of plane waves), 
while the local mean height itself varies slowly on the scale of the mean wavelength 
and mean period. This approximation is well justified, since the envelope A{X, Y, T) in 
equation (l23ll is slowly varying for Ak/k^ <^ 1 and <^ 1, while the higher harmonics 
A2{X, Y, T) and A2,{X, Y, T) are suppressed by factors of e and e^, respectively. Taking 
the local mean intensity to be distributed, and convolving the distribution of the 
mean intensity with the Rayleigh distribution around the mean intensity, we obtain as 
in the linear case a K-distribution (fT3l) for the total distribution of wave heights. 

In figure [HI, each data set is fit to the K-distribution of equation f[T3l) . arising from 
the local Rayleigh approximation. We see that the fits, indicated by dashed and dotted 
lines, perform adequately for probabilities down to 10~®, where statistical noise begins to 
dominate. In particular, we clearly observe the crossover between the Gaussian behavior 
([Hj) at small to moderate heights and the asymptotic exponential behavior ( JT71) at large 
heights. However, systematic deviations do exist, which are especially visible at larger 
values of e, corresponding to smaller values of the N (degrees of freedom) parameter. 
These systematic deviations are in large part due to the fact that the true wave height 
distribution for any given set of input parameters exhibits spatial dependence, evolving 
from the original Rayleigh distribution imposed by incoming boundary conditions to 
the broader K-distribution, and then gradually back to a Rayleigh distribution as the 
wave energy is transferred to longer wavelengths and the steepness decreases [42]. An 
example of this spatial dependence appears below in figure [Til Thus, a more accurate 
model for the total wave height distribution consists of a sum of several K-distributions, 
or equivalently the tail of the full distribution may be modeled by a K-distribution 
multiplied by a prefactor C < 1, as discussed in reference [13]. Nevertheless, as seen in 
figure [HI equation ( [T3l) correctly describes wave height probabilities at the ±20% level 
of accuracy, allows for an extremely simple one-parameter characterization of the wave 
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Figure 8. The distribution of wave heights, in units of the significant wave height, is 
calculated for three nonzero values of the steepness e (upper three solid curves), and 
compared with the random seas model of equation (lowest solid curve). In each 
case, the dashed or dotted curve is a best fit to the K-distribution of equation 
Here the we fix the angular spread — 2.6° and wave number spread Afc/feg = 0.1 
of the incoming sea. 



height distribution, and facihtates easy comparison between the effects of hnear and 
nonhnear focusing. 

3.3. Scaling with Input Parameters 

Given the single-parameter approximation of equation ( IT3l) . it is sufficient to explore 
the dependence of the parameter N on the input variables describing the incoming sea, 
specifically the initial angular spread A^, the initial wave number spread Ak/ko, and 
the initial steepness e. In the two panels of figure [H we fix the steepness at 5 = 0.032 
and show the scaling of with A6 and Sk/ko, respectively. Given that the Benjamin- 
Feir instability for a monochromatic wave in one dimension [H] is at the root of the 
nonlinear instability in the general case, it is not surprising that stronger deviations 
from the Rayleigh model, as indicated by smaller values of A^, occur as A6 or Ak is 
reduced, consistent with earlier results [l5l HGj HTJ HH] . Specifically, we find 



where a, 6 ~ 1. 

We note that the scaling of with incoming angular spread A6 for nonlinear 
focusing, ~ A^, is only half as strong as the scaling A^ ~ (A0)^ arising from 
linear wave scattering by currents, as implied by equations ([6]) and ([8]). Thus, smaller 
angular spreads A^ are needed for the nonlinear focusing mechanism to be effective, as 
compared with linear focusing by currents. This is easily seen by comparing the range 
of A^ in figure [9] with the corresponding range in table [1] for the linear mechanism. On 




(27) 
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Figure 9. The best-fit N value (equation describing tlie wave height probability 
distribution is shown as a function of the initial angular spread A6 and initial wave 
number spread Ak/ko of the incoming sea. The steepness is fixed at e = 0.032. The 
left panel shows the scaling of N with A9, with the line showing the best-fit scaling 
N ^ (A6')i °* for Ak/ko = 0.15. The right panel shows the scaling of N with Ak/ko, 
with the line showing the best-fit scaling N ^ (Ak/ko)^-^^ for A9 — 2.6°. 



the other hand, figure [9] and equation (!27j) both imply that the nonhnear mechanism 
exhibits significant sensitivity to the spectral width Ak/k^, consistent with previous 
findings [6l HHl [501 EU [52], [53] . This is to be contrasted with the linear mechanism of 
rogue wave formation, which is insensitive to the spectral width at leading order in 
Ak/ko M- 




1 -_ 
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Figure 10. The best-fit N value (equation (|13p ) describing the wave height probability 
distribution is shown as a function of the steepness e for several values of the initial 
angular spread A9. The initial wave number spread is fixed at Ak/ko = 0-lj as in 
figure[Hl The line shows the best-fit scahng N ~ e''^-^ for A0 = 2.6°. 

Finally, in figure [TOl we fix the steepness Ak/ko = 0.1, as in figure [HI and examine 
the scaling of the value with the steepness e, for several values of the initial angular 
spread A^. As e grows, N decreases, indicating greater deviations from the Rayleigh 
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distribution. Again, we observe good power-law scaling with the steepness in the range 
of parameters considered here. We have 

N ^e" (28) 

where c ~ —3. At larger values of the steepness (not shown), saturation occurs. 



A^ 


^(2.2) 


^(3.0) 


2 


1.1 • 10^ 


5.2 ■ 10^ 


5 


37 


7.3 ■ 10^ 


10 


16 


1.3 ■ 103 


20 


6.8 


2.2 • 10^ 


50 


2.9 


27 


100 


1.8 


7.8 



Table 2. The enhancement in the probabihty of rogue wave formation (wave height 
2H — 2.2 SWH) as well as the enhancement of the probability of extreme rogue wave 
formation (wave height 2H — 3.0 SWH) are calculated for several values of the N 
parameter, as in table [TJ 

Table [21 calculated analogously to tabled] in the previous section, aids in extracting 
the implications of figures |9] and [10] by indicating the quantitative relationship between 
the A^ value and the enhancement in rogue wave and extreme rogue wave occurrence 
probabilities. We note that even at A^ values between 50 and 100, corresponding to 
the upper range of values in figures [9] and [THl the occurrence of extreme rogue waves is 
enhanced by an order of magnitude. Exponentially larger enhancement is predicted for 
parameters associated with smaller values of A^. 

4. Combined Effect of Nonlinear and Linear Focusing 

Finally, we discuss the possibility of even greater enhancement in the rogue wave 
formation probability when linear and nonlinear mechanisms are acting together |l2l |13] . 
In this context, it is important to consider again the spatial scales associated with rogue 
wave development in the two mechanisms. We recall that when an incoming random 
sea is linearly scattered by strong currents, the first singularities in the ray dynamics 
occur after a distance scale L ~ ^{uj:ms/v)~'^^'^, as discussed in section 12.11 These 
first singularities are the ones associated with the highest probability of rogue wave 
formation, as subsequent random scattering exponentially stretches the phase space 
manifold and leads to ever smaller density associated with each second- and higher-order 
singularity [iHl [H]. At distances S> L, the pattern of hot and cold spots becomes less 
and less prominent, the ray density again becomes nearly uniform (see figure [T|), and the 
wave height distribution asymptotically approaches again the Rayleigh limit. Similarly, 
nonlinear evolution as described by the NLSE without current (equation (^1^ with 
$c = 0) occurs on a typical distance scale 1/ke^. On distance scales larger than 1/ke^, 
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energy transfer from smaller to larger wavelengths (i.e., the frequency downshifting effect 
mentioned previously in Sec. 13. ip results eventually in a decline in the steepness and 
again an approach towards the limiting Rayleigh distribution [M| [T0| [T2] . 
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Figure 11. The fourth moment of the wave height distribution is shown as a function 
of evolution distance, starting in each case from a Longuet-Higgins random sea with 
mean wave speed v — 7.81 m/s, initial angular spread — 5.2°, and wave number 
spread tJ^jk^ = 0.1. The three situations considered are: (a) linear scattering by 
random currents of rms speed itj-ms = 0.2 m/s and eddy correlation length — 800 m, 
(b) nonlinear evolution with initial steepness e — 0.032 and without currents, and (c) 
nonlinear evolution in the presence of currents. 



This behavior is illustrated in figure [TT] for linear evolution with random currents, 
nonlinear evolution in the absence of currents, and for a scenario in which the two 
mechanisms are both active. Here we use the fourth moment H'^ as a convenient 
measure of the size of the tail of the wave height distribution. Note that for the 
chosen parameters, the distance scales associated with linear and nonlinear rogue wave 
formation are comparable. Clearly, in this case currents have a greater effect than 
nonlinear focusing, but the strongest deviations from Rayleigh statistics are observed 
when linear scattering and nonlinear interaction are both present. 

The total wave height distributions for these same three scenarios, and the 
probability enhancement over the predictions of the Longuet-Higgins model, are shown 
in figure [121 As noted above in section 13.21 when wave height data is collected over a 
large spatial field that includes some areas of very strong deviations of Rayleigh statistics 
and other areas where such deviations have not yet had an opportunity to develop, the 
full distribution may not be well approximated by a single K-distribution, but the tail 
may still be well approximated in this way, since it is dominated by data from those 
areas where deviations are strongest ^43j. This is indeed what we clearly observe in 
figure [T2| for the scenario where nonlinearity and currents are both present. 

Again we see from figure [12] that deviations from Rayleigh statistics become ever 
more pronounced as taller and taller waves are considered, as expected from the 
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Figure 12. Left panel: The cumulative distribution of wave heights, in units of the 
significant wave height, is obtained for the same three scenarios as are considered in 
figure [11] In each case, the solid curve is a fit to a K-distribution with (a) = 16 
for linear scattering by currents, (b) = 29 for nonlinear evolution, and (c) N — 5.1 
when linear and nonlinear focusing are acting in concert. The Rayleigh distribution 
(A^ = oo) is shown for reference. Right panel: In each of the three scenarios, the 
probability enhancement factor -Ptotai(-ff)/^'Rayicigh(-ff) is obtained from the data. 

asymptotic form of the K-distribution (equation ( !T7|) ). In particular, in this example we 
see that the probability of forming an extreme rogue wave (wave height = 3 SWH) is 
enhanced by a factor 90 due to nonlinear interaction, by a factor of 380 due to focusing 
by currents, and by a factor of 2600 when the two mechanisms are combined. 

5. Conclusions and Outlook 

It will take some time to sort out the mechanisms of rogue wave formation with complete 
certainty. All potentially important factors and mechanisms ought to be included in the 
discourse, which we hope will someday lead to agreement about the several formation 
mechanisms and their interactions. More importantly, predictive tools leading to safer 
navigation should eventually emerge. One of the seemingly important factors, which 
might be called "statistical focusing," is highlighted here. In terms of wave propagation, 
statistical focusing is a linear effect (although it is nonlinear dynamics at the level of 
ray tracing). It leads to large enhancements in the frequency of rogue wave formation 
under reasonable sea state assumptions. 

Statistical focusing combines the effects of deterministic wave refraction by current 
eddies with Longuet-Higgins statistical ideas under realistic conditions. The key notion 
is that the focusing effects of eddies, which would be very dramatic on an (unrealistic) 
monochromatic and unidirectional not altogether washed out when realistic 

frequency and directional dispersion are included. Essentially, deterministic caustics 
present in the unrealistic idealization are smoothed into hot spots, which are then treated 
statistically within Longuet-Higgins theory. The hot spots dominate the statistics 
in the tail of the wave height distribution. This amounts to a nonuniform sampling 
version of Longuet-Higgins theory, with a solid basis for the nonuniform energy density 
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distributions used. 

Since nonlinear effects are also important, we have examined them alone within 
the popular fourth-order nonlinear Schrodinger equation (NLSE) approximation for 
nonlinear wave evolution under realistic seaway conditions. Finally, we have investigated 
the combined effect of nonlinear wave evolution and statistical focusing. We find 
that strongest deviations from Rayleigh statistics are observed when linear scattering 
(statistical focusing) and nonlinear interaction (NLSE) are both present. However, 
for the parameters chosen here at least, the linear scattering due to eddies was more 
important than the nonlinear effects, which require large steepness or a very narrow 
range of propagation directions to become significant. 

We have presented a measure closely related to the probability of rogue wave 
formation, the freak index 7. This could conceivably become the basis for a probabilistic 
forecast of rogue wave formation, in the spirit of rainfall forecasts. 

There are at least three clear directions for future development of the work presented 
here. First, both the computer simulations and the theory must be developed further to 
explore fully and systematically the combined effects of nonlinear and linear focusing. 
This will also involve investigating in depth the underlying mechanism through which 
the formation of hot and cold spots is aided by nonlinear focusing. Secondly, a better 
understanding is needed of the stability of the hot spot patterns under slow changes in 
the current field or in the spectrum or directionality of the incoming sea. The strength 
of what might be called scintillation or twinkling [25] in analogy with the case of light 
traveling through the atmosphere will have important consequences for the predictive 
power of the model. Thirdly, and most importantly, there is a clear need to compare the 
model simulations with observations and experiments. Although comprehensive global 
data are not available at this point, it may be possible to compare the results to local 
observations where data are more readily available, e.g., in the North Sea. 

Whatever the final word is on rogue wave formation (or final words, because there 
may be more than one mechanism), it must involve a reallocation of energy from a larger 
area to a smaller one. Waves cannot propagate and increase in height at no expense 
to their neighbors: the energy has to come from somewhere, and the effect must be to 
reduce the wave energy somewhere else. The focusing mechanism is clear in this respect: 
hot spots form and cold spots do too, according to a ray tracing analysis, maintaining 
energy balance |55j . 
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